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<^ | Abstract 
(N 

We describe a dynamic programming algorithm for exact counting and exact uni- 
form sampling of matrices with specified row and column sums. The algorithm runs 
in polynomial time when the column sums are bounded. Binary or non-negative inte- 
ger matrices are handled. The method is distinguished by applicability to non-regular 
margins, tractability on large matrices, and the capacity for exact sampling. 
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: 1 Introduction 

q 

^ \ Let iV(p, q) be the number of m x n binary matrices with margins (row and column sums) 

p = (px, . . . ,p m ) G N m , q = (qi, . . . , q n ) G N n respectively, and let M(p, q) be the 
corresponding number of N- valued matrices. In this paper we develop a technique for 
efficiently finding iV(p, q) and M(p, q). Uniform sampling from these sets of matrices is 
^ ! an important problem in statistics 0, and the method given here permits efficient exact 

uniform sampling once the underlying enumeration problem has been solved. 

Since a bipartite graph with degree sequences p = . . . , p m ) G N"\ q = (gi, . . . , q n ) G 
N n (and m, n vertices in each part respectively) can be viewed as a m x n matrix with 
row and column sums (p, q), our technique applies equally well to counting and uniformly 
sampling such bipartite graphs. Under this correspondence, simple graphs correspond to 
binary matrices, and multi graphs correspond to N-valued matrices. 
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The distinguishing characteristic of the method is its tractability on matrices of non-trivial 
size. In general, computing M(p, q) is #P-complete [10], and perhaps N(p, q) is as well. 
However, if we assume a bound on the column sums then our algorithm computes both 
numbers in polynomial time. After enumeration, uniform samples may be drawn in poly- 
nomial expected time for bounded column sums. To our knowledge, all previous algorithms 
for the non-regular case require super-polynomial time (in the worst case) to compute these 
numbers, even for bounded column sums. (We assume a description length of at least m+n 
and no more than mloga + nlogb, where a = maxpi, b = maxgj.) In general (without 
assuming a bound on the column sums), our algorithm computes N(p, q) or M(p, q) in 
0(m(ab + c)(a + 6) 6_1 (6 + c) 6_1 (logc) 3 ) time for m x n matrices, where a = m&xpi, 
b = max qi, and c = J^Pi = Yl Qi- After enumeration, uniform samples may be drawn in 
O ( mc log c) expected time. 

In complement to most approaches to computing M(p, q), which are efficient for small 
matrices with large margins, our algorithm is efficient for large matrices with small 
margins. For instance, in Section |4] we count the 100 x 100 matrices with margins 

(70, 30, 20, 10, 5( 6 \ 4(10) ; 3(20) ^ 2 (60)) 5 (4(80^ 3(20)) (where x {n) denotes x repeated n times). 

To illustrate the problem at hand, consider a trivial example: if p = (2,2, 1, 1), q = 
(3, 2, 1), then iV(p, q) = 8 and M(p, q) = 24. The 8 binary matrices are below. 
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The paper will proceed as follows: 

§2 Main results 

§3 Brief review 

§4 Applications 

§5 Proof of recursions 

§6 Proof of bounds on computation time. 



2 Main results: Recursions, Bounds, Algorithms 

Introducing the following notation will be useful. Taking N := {0, 1,2,...}, we consider 
N n to be the subset of N°° := {(r u r 2 , . . . ) : n e N for i = 1, 2, . . . } such that all but 
the first n components are zero. Let L : N°° — » N°° denote the left- shift map: Lr = 
(r 2 , . . . , r n , 0, 0, . . . ). Given r, s e N™, let r\s := r — s + Ls, (which may be read as "r 
reduce s"), let 
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and let f denote the vector of counts, f := (f i, f 2 , . . . ) where := #{j : rj = i}. We write 
r < s if Ti < s,i for all i. Given n G N, let C n (k) := {r 6 N" : £V r*j = k} be the n-part 
compositions (including zero) of k, and given s G N", let C s (k) := {r G C n (k) : r < s}. 
For (p, q) G N m x N n , define the numbers 



N{p, q) := #{X G {0, l} mxn : J^xy = p*, J^xy = fc, for 1 < i < m, 1 < j < n}, 



M(p, q) := #{X G N mx " : £ x - = Pij ^ x y = g„ for 1 < i < m, 1 < j < n}. 



Since N(p, q) and M(p, q) are fixed under permutations of the row sums p and column 
sums q, and since zero margins do not affect the number of matrices and can effectively 
be ignored, then we may define N(p, q) := N(p, q) and M(p,q) := M(p, q) without 
ambiguity. We can now state our main results. 

Theorem 2.1 (Recursions) The number of matrices with margins (p, q) G W n x N n is 

given by 



s€C r + is (pi) V 7 

where r = q, and in (2), we sum over all s such that s G C r+Ls (pi). 

Proofs will be given in Section[5] The Gale-Ryser conditions lfTTll29l simplify computation 
of the sum in (1) by providing a necessary and sufficient condition for there to exist a binary 
matrix with margins (p, q): if q[ := #{j : qj > i} and pi > • • • > p m , then N(p, q) ^ 
if and only if Yli=iPi — Si=i9i for all j < m and ^iliPi = Y^i=iQi- This is easily 
translated into a similar condition in terms of (p, q) and N(p, q). The following recursive 
procedure can be used to compute either N(p, q) or M(p, q). 

Algorithm 2.2 (Enumeration) 

Input: (p, q), where (p, q) G N m x N n are row and column sums such that J2i Pi = J2i Qi- 

Output: N(p, q) (or M(p, q) ), the number of binary (or N-valued) matrices. 

Storage: Lookup table of cached results, initialized with N(0, 0) = 1 (or M(0, 0) = 1). 

(1) IfN(p, q) is in the lookup table, return the result. 

(2) In the binary case, if Gale-Ryser gives N(p, q) = 0, cache the result and return 0. 

(3) Evaluate the sum in Theorem \2.1\ recursing to step (I) for each term. 

(4) Cache the result and return it. 



j 




for binary matrices, and 




for N-valued matrices, 
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Let T(p,q) be the time (number of machine operations) required by Algorithm 12.21 to 
compute N(p, q) or M(p, q), after performing an 0(n 3 ) preprocessing step to compute all 
needed binomial coefficients. (It turns out that computing M(p, q) always takes longer, but 
the bounds we prove apply to both N(p, q) and M(p, q).) We give a series of bounds on 
T(p, q) ranging from tighter but more complicated, to more crude but simpler. The bounds 
will absorb the 0(n 3 ) pre-computation except for the trivial case when the maximum col- 
umn sum is 1. 



Theorem 2.3 (Bounds) Suppose (p, q) G N m x N n , a = m&xpi, b = maxg^, and c = 
Y.Pi = Yl ft- Then 

(1) r(p, q) < 0((«6 + c)(log cf ± (« l h _ ~ X ) ^ + • • • + »J + h ~y 

(2) T(p,q) < 0(m(a6 + c)(a + 6) 6 - 1 (6 + c)^ 1 (logc) 3 ), 

(3) T(p,q) < 0{mn 2b - l (\ognf) for bounded b, 

(4) T(p, q) < 0(mn b (log n) 3 ) for bounded a, b. 



Remark Since we may swap the row sums with the column sums without changing 
the number of matrices, we could use Algorithm 12.21 on (q, p) to compute iV(p, q) or 
M(p, q) using T(q, p) operations, which, for example, is 0(nm a (logm) 3 ) for bounded 
a, b. T(p, q) also depends on the ordering of the row sums pi, . . . ,p m as suggested by 
Theorem 12.3( 1). and we find that putting them in decreasing order p 1 > ■ ■ ■ > p m tends 
to work well. Algorithm 12.21 is typically made significantly more efficient by using the 
Gale-Ryser conditions, and this is not accounted for in these bounds. Although we observe 
empirically that this reduces computation tremendously, we do not have a proof of this. 



Algorithm I2.2l traverses a directed acyclic graph in which each node represents a distinct set 
of input arguments to the algorithm, such as (p, q). Node (u, v) is the child of node (p, q) 
if the algorithm is called (recursively) with arguments (u, v) while executing a call with 
arguments (p, q). If the initial input arguments are (p, q), then all nodes are descendents of 
node (p, q). Meanwhile, all nodes are ancestors of node (0, 0). Note the correspondence 
between the children of a node (u, v) and the compositions s G C v (u\) in the binary 
case, and s G C v+Ls (wi) in the N- valued case, under which s corresponds with the child 
(Lu, v\s). We also associate with each node its count, the number of matrices with the 
corresponding margins. 

As an additional benefit of caching the counts in a lookup table (as in Algorithm 12.21 ), once 
the enumeration is complete we obtain an efficient algorithm for uniform sampling from 
the set of (p, q) matrices (binary or N-valued). It is straightforward to prove that since the 
counts are exact, the following algorithm yields a sample from the uniform distribution. 
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Algorithm 2.4 (Sampling) 

Input: 

■ Row and column sums (p, q) G N m x N n such that YliPi = J2i Qi- 

■ Lookup table of counts generated by Algorithm \2.2\ on input (p, q). 

Output: A binary (or N-valued) matrix with margins (p, q), drawn uniformly at random. 

(1) Initialize (u, v) <— (p, q). 

(2) If{u,v) = (0,0), exit. 

(3) Choose a child (Lu, v\s) 0/(11, v) wzY/i probability proportional to its count times the 
number of corresponding rows (that is, the rows r G C v (ui) such thatv — r = v\s.J 

(4) Choose a row uniformly among the corresponding rows. 

(5) (u,v) <- (Lu,v- r). 

(6) Goto (2). 

In step ©, there are ( v ) corresponding rows r in the binary case, and (*+ Ls ) in the N- 
valued case. In step (0), in the binary case of course we only choose among r G {0, l} n . In 
Section [6] we prove that Algorithm I2.4l takes O(mclogc) expected time per sample, where 



3 Brief review 

We briefly cover the previous work on this problem. This review is not exhaustive, focusing 
instead on those results which are particularly significant or closely related to the present 
work. Let H n (r) and H*(r) denote M(p, q) and N(p, q), respectively, when p = q = 
(r, . . . , r) G N n . The predominant focus has been on the regular cases H n (r) and H*(r). 

Work on counting these matrices goes back at least as far as MacMahon, who applied 
his expansive theory to find the polynomial for H 3 (r) [|2T| (Vol II, p. 161), and developed 
the theory of Hammond operators, which we will use below. Redfield's theorem [|2~8l . 
inspired by MacMahon, can be used to derive summations for some special cases, such as 
H n (r), H*(r) for r = 2, 3, and in similar work Read [|26l l27l used Polya theory to derive 
these summations for r = 3. Two beautiful theoretical results must also be mentioned: 
Stanley OTI proved that for fixed n, H n (r) is a polynomial in r, and Gessel lfT2~l[T3l showed 
that for fixed r, both H n (r) and H*(r) are P-recursive in n, vastly generalizing the linear 
recursions for H n (2), H*(2) found by Anand, Dumir, and Gupta JTJ. 

We turn next to algorithmic results more closely related to the present work. McKay 
11221 121 has demonstrated a coefficient extraction technique for computing N(p, q) in the 
semi-regular case (in which p = (a, . . . , a) G N m and q = (6, . . . , b) G W 1 ). To our 
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knowledge, McKay's is the most efficient method known previously for N(p, q). By our 
analysis it requires at least VL{mn b ) time for bounded a, b, while the method presented here 
is 0(mn b (\ogn) 3 ) in this case. Since this latter bound is quite crude, we expect that our 
method should have comparable or better performance, and indeed empirically we find that 
typically it is more efficient. If only b is bounded, McKay's algorithm is still VL{mn b ), 
but the bound on our performance increases to 0(mn 26_1 (logn) 3 ), so it is possible that 
McKay's algorithm will outperform ours in these cases. Nonetheless, it is important to 
bear in mind that McKay's algorithm is efficient only in the semi-regular case (while our 
method permits non-regular margins). If neither a nor b is bounded, McKay's method is 
exponential in b (as is ours). 

Regarding M(p, q), one of the most efficient algorithms known to date is LattE (Lattice 
point Enumeration) |fT9l , which uses Barvinok's algorithm [J21 to count lattice points con- 
tained in convex polyhedra. It runs in polynomial time for any fixed dimension, and as a 
result it can compute M(p, q) for astoundingly large margins, provided that m and n are 
small. However, since the computation time grows very quickly with the dimension, LattE 
is currently inapplicable when m and n are larger than 6. There are similar algorithms 
Il23ll20l l3l that are efficient for small matrices. 

In addition, several other algorithms have been presented for finding iV(p, q) (such as 
lfT8l [321 [331 [24ll ) and M(p, q) (see review [[El) allowing non-regular margins, however, 
it appears that all are exponential in the size of the matrix, even for bounded margins. 
While in this work we are concerned solely with exact results, we note that many useful 
approximations for iV(p, q) and M(p, q) (in the general case) have been found, as well as 
approximate sampling algorithms [[T7l l7l[l4ll4l[l6l. 

4 Applications 

4.1 Occurrence matrices from ecology 

The need to count and sample occurrence matrices (binary matrices indicating observed 
pairings of elements of two sets) arises in ecology. A standard dataset of this type is "Dar- 
win's finch data", a 13 x 17 matrix indicating which of 13 species of finches inhabit which of 
17 of the Galapagos Islands. The margins of this matrix are (14, 13, 14, 10, 12, 2, 10, 1, 10, 
11,6,2, 17), (4, 4, 11, 10, 10,8,9, 10,8,9,3, 10, 4,7,9,3,3). We count the number of such 
matrices to be 67,149,106,137,567,626 (in 1.5 seconds) confirming 0. Further, we sample 
exactly from the uniform distribution over this set at a rate of 0.001 seconds per sample. 
(All computations were performed on a 64-bit 2.8 GHz machine with 6 GB of RAM.) A 
similar dataset describes the distribution of 23 land birds on the 15 southern islands in the 
Gulf of California [6): this binary matrix has margins (14, 14, 14, 12,5, 13,9, 11, 11, 11, 
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1 1, 1 1, 7, 8, 8, 7, 2, 4, 2, 3, 2, 2, 2), (21, 19, 18, 19, 14, 15, 12, 15, 12, 12, 12, 5, 4, 4, 1), for 
which we count 839,926,782,939,601,640 corresponding binary matrices. Counting takes 1 
second, and sampling is 0.002 seconds per sample. One more example of this type: for bird 
species on the California Islands fl25l we find that there are 1 ,360,641,571,195,21 1,109,388 
binary matrices with margins (1, 4, 3, 2, 1, 1, 1,5, 1,3, 1, 4, 4, 5, 1,2, 1, 5, 4, 5, 3, 7, 1, 
3, 2, 4, 1, 3, 2, 4, 6), (2, 14, 24, 8, 2, 5, 20, 15) in 4 seconds; samples take 0.003 seconds 
each. Larger matrices can be handled as well, provided the margins are small. For exam- 
ple, we count 860585058801817078819959949756.. .000 (459 digits total, see Appendix) 
100 x 100 matrices with margins (70, 30, 20, 10, 5< 6 \ 4( 10 ), 3( 20) , 2( 60 )), (4( 80 \ 3< 20 >) (where 
x^ denotes x repeated n times) in 46 minutes. We know of no previous algorithm capable 
of efficiently and exactly counting and sampling from sets such as this. 

4.2 Ehrhart polynomials of the Birkhoff polytope 

Stanley OTI proved a remarkable conjecture of Anand, Dumir, and Gupta [1J: given n G N, 
H n {r) is a polynomial in r (where H n {r) = M(p, q) with p = q = (r, r, . . . , r) G N n ). 
Given H n (l), . . . , iJ n (( n ~ 1 )), one can solve for the coefficients of H n (r) (as we describe 
below). These polynomials have been computed for n < 9 by Beck and Pixton (3). As an 
application of our method, we computed them for n < 8, and found that the computation 
time is comparable to that of Beck and Pixton. For n = 4, . . . , 8, the numbers H n (r) for 
r = 1, . . . , ( n ~ are listed in the Appendix, and the polynomial H±(r) is displayed here as 
an example. Our results confirm those of Beck and Pixton. 

H 4 (r) = 1 + (65/18)r + (379/63)r 2 + (35117/5670)r 3 + (43/10)r 4 
+ (1109/540)r 5 + (2/3)r 6 + (19/135)r 7 + (ll/630)r 8 + (ll/11340)r 9 . 

The coefficients of H n {r) can be determined by the following method. By Stanley's 
theorem [31], H n {r) is a polynomial in r such that (a) deg H n (r) = (n — l) 2 , (b) 
H n (-1) = ■■■ = H n (-n+l) = 0,and(c) H n {-n-r) = (-l)^" 1 ) H n (r) for r G N. For 
each n — 4, . . . , 8, we perform the following computation. Let k = ( n ~ x ) and d = (n— l) 2 . 
Compute the numbers H n (r) for r = 0,1, ... ,k using Algorithm 12.21 and form the vec- 
tor v := (H n (-n - k + 1), . . . , H n (k)) T G Z d+1 using (b) and (c). Form the matrix 
A = ((i - n - ky- 1 )^ G z( d+1 ) x ( d+1 ), and compute u = A _1 v. Then by (a), 




3=0 
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4.3 Contingency Tables 

As an example of counting contingency tables with non-regular margins, we count 
620017488391049592297896956531. ..000 (483 digits total, see Appendix) 100 x 100 ma- 
trices with margins (70, 30, 20, 10, 3^ 20 \ 2^), 3^) (where x^ denotes 
x repeated n times) in 118 minutes. Again, we know of no previous algorithm capable 
of efficiently and exactly counting and sampling from sets such as this. (However, for 
small contingency tables with large margins, our algorithm is much less efficient than other 
methods such as LattE.) Exact uniform sampling is possible for contingency tables as well, 
which occasionally finds use in statistics |[9]|. 

5 Proof of recursions 

We give two proofs of Theorem 12.11 The first is a "direct" proof, which provides the basis 
for the sampling algorithm outlined above. In addition to the direct proof, we also provide a 
proof using generating functions which is seen to be a natural consequence of MacMahon's 
development EH of symmetric functions, and yields results of a more general nature. 

5.1 Preliminary observations 

For r G N n , let r' denote the conjugate of r, that is, r[ = #{j : r,j > i} for i = 1, 2, 3, 

For r, s G N°°, let r A s denote the component- wise minimum, that is, (ri Asi, r-i AS2, ■ ■ ■ ). 
In particular, r A 1 = (r\ A 1, r 2 A 1, . . . ). Recall our convention that N n is considered to be 
the subset of N°° such that all but the first n components are zero. (Similarly, we consider 

Z n C Z°°.) 

Lemma 5.1 Let u,veN" such that u < v. 

(1) Suppose s G N d for some d G N. Then v — u = v\s if and only ifs = v' — (v — u)'. 

(2) If s — v' — (v — u)' then ^2 Si = Yl u i an d s < v + Ls. 

(3) Ifs = v' — (v — u)' and u < v A 1 then s < v. 

Proof (1) Letting I be the identity operator, a straightforward calculation shows that for 
any d G N, r G Z d , we have (J2T=o Lk ) ( J - L)t = t and (I - L) (J2T=o L>c ) r = r > that 
is, (I — L)^ 1 = J2T=o L k on where I is the identity operator. Further, (I — L)~ l r = r'. 
Thus, v — u = v\s = v — s + Ls if and only if (I — L)s = v — v — u if and only if 
s = (I - Ly\v - V=u) = v' - (v - u)'. 
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(2) Ifs = v'-(v-u)'thenX)si = XX'-XX V_U X = X/^XX^-^) = E u i since 
Yl r 'i = Yl r i f° r a ^ r G N n . By (1), v — s + Ls = v\s = v — u > and so s < v + Ls. 

(3) If s = v' — (v — u)' and u < v A 1 then by definition, Sj = #{j : Vj > i} — #{j : 
Vj - uj >i} = #{j : Vj = i and Uj = 1} < #{j : Vj = i] = fy. | 

Lemma 5.2 Le? v G N n , fceN, and to /(u) = v' - (v - u)' for u G N" such that u < v. 

(1) /(C vA1 (A;)) = C*{k), and for any s G C*(k), #{u G C vA1 (A;) : /(u) = s} = (*). 

(2) f{C v (k)) = {s : s G C^ +is (£;)}, and for any s swc/z ffczf s G C*+ is (A;), #{u G 



Proof (1) /(C vA1 (A;)) c C^(Jfe) follows from Lemma |5J22 and 3). Let s G C^(ife). 
Choose u as follows. For i = 1, 2, 3, ... , choose Sj of the positions j such that Vj = i, 
and set Uj = 1 for each chosen j. (Set Uj = for all remaining j.) This determines some 
u G C vA1 (/c) such that Sj = #{j : Vj — i and itj : = 1} for all z. Furthermore, it is not hard 
to see that any such u is obtained by such a sequence of choices. Now, as in the proof of 
Lemma I5TTL 3), = #{j : Vj = i and u, } ■ = 1} if and only if /(u) = s (when u < v A 1). 
Hence, f(C vA1 (k)) D C v (k), and since there were ( v ) possible ways to choose u, then 
this proves (1). 

(2) f(C v (k)) C {s : s G C* +Ls (k)} follows from Lemma EJJ2). Suppose s G C* +Ls (k). 
Let r = v and t = r\s. Note that t > since s < r + Ls. Also, r, s, t G N d where 
d = m&xVj. Choose u as follows. First, consider the positions j in v at which Vj = d. 
There are Q d ) ways to choose td of these positions. Having made such a choice, we 
set Uj = Vj — d = for each such j that was chosen. Next, consider the r^^i positions 
j at which Vj = d — 1, in addition to the — ^ remaining positions at which Vj = d. 
There are ( rd ~ 1 ^ rd ~ <d ' ) ) ways to choose td-i of these. Having made such a choice, we 
set Uj = Vj — (d — 1) for each such j that was chosen. Continuing in this way, for i = 
d — 2, . . . , 1: consider the positions j in v which Vj = i, in addition to the r i+1 + • • • + 
r d — td — ■ ■ ■ — U+i remaining positions at which vj > i, choose t { of these (in one of 



was chosen. After following these steps for each i, set Uj = Vj for any remaining positions 
j. This determines some u such that < u < v. 

Now, for i — d, d — 1, . . . , 1, we have chosen ti positions j and we have set Uj = Vj — i. 
That is, ti = #{j : Vj — Uj = i}, and so t = v — u. Hence, v — u = v\s (by the definition 



of t), so s = /(u) by Lemma l5TTf 1). and additionally, X u j = X s j = k by 15. 1 T 2V Thus, 
we have shown that f(C v (k)) D {s : s G C*+ Ls (k)}. ' 



C{k) : /(u) = s} = (^ s ). 




ways), and set Uj 



Vj — i for each such j that 
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Using tj = r,j — Sj + Sj + \ (the definition of t), we see that there were 

V<A / r d _i + (r d - t d )\ fr 1 + r 2 J i V r d - t d - ■ 



tdj \ td-i J \ ti 

r d \ ( Td-\ + s d \ (r 1 + s 2 \ fr + Ls 

Sd) V S d~l J V s ! 



> 



ways to make such a sequence of choices, where the inequality holds since s < r + Ls. 
Hence, there are at least ( r+ s Ls ) distinct choices of u G C v (k) such that /(u) = s. On the 
other hand, given any u G C v (k) such that /(u) = s, we have t = v — u (by Lemma 
I5.HT )). thus ti = : Uj = Vj — %}, and since Vj > i for any j such that Uj = Vj — i, such 
a u is obtained by one of the sequences of choices above. Hence, #{u G C v (k) : /(u) = 

s} = r s Ls ). i 



5.2 Direct proof 

We are now prepared to prove Theorem 12.11 Recall the statement of the theorem: 
The number of matrices with margins (p, q) G N m x N n is given by 

(1) N(p,r) = I jiV(Lp, r\s) for binary matrices, and 

seC r (pi) ^ 

M(Lp,r\s) for N-valued matrices, 

where r = q, and in (2), we sum over all s such that s G C r+Ls (pi). 
Proof of Theorem 2.1 

(1) First, we prove the binary case. Let (p, q) G N m x N n , r = q. Using Lemma [5T2T 1), 
define the surjection / : C qA1 (pi) — > C r (pi) by /(u) = q' — (q — u)'. Then 

iV(p,r) = iV(p,q)i N{Lp,q-n) 

= E E ^p.q-u)= E (*Wp>*vo- 

Step (a) follows from partitioning the set of (p, q) matrices according to the first row u G 
C qA1 (Pi) of the matrix. Step (b) partitions C qA1 (pi) into the level sets of /, that is, the sets 
/-i(s) = {u G C qA1 (pi) : /(u) = s} as s ranges over /(C qA1 (pi)) = C r ( Pl ). Step (c) 
follows since if /(u) = s then q — u = r\s (by Lemma I51T 1)) and thus N(Lp, q — u) = 
N(Lp, r\s), and since #/" 1 (s) = (*) (by LemmaQl)) • This proves 0;i). 
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(2) Now, we consider the N-valued case. Let S = {s : s G C r+Ls (pi)}. Using Lemma 
I5.2f 2), define the surjection g : C q (pi) — > S by g(u) = q' — (q — u)'. Then, similarly, 

M(p,r) = M(p,q) ( = ) ^ M(Lp, q - u) 

ueci(pi) 

£ M(Lp,q-u)iW r+ s iS Wp,r\s). 

As before, step (a) follows from partitioning the set of matrices according to the first row 
u G C q (pi), step (b) partitions C q (pi) into the level sets of g, and step (c) follows since 
^^(s) = ( r + Ls ) (by LemmaE22)). This proves TheoremO | 

5.3 Generating function proof 

In addition to the direct approach above, one may also view the recursions as the application 
of a certain differential operator to a certain symmetric functions. Although such operators 
were used extensively by MacMahon ETTl on problems of this type, at first it would appear 
that for computation this approach would be hopelessly inefficient in all but the simplest 
examples. In fact, it turns out that a simple observation allows one to exploit regularities 
in the present problem, reducing the computation time to polynomial for bounded margins. 
Specifically, when there are many columns with the same sum, the symmetric function un- 
der consideration has many repeated factors, and the action of the operator in this situation 
takes a simplified form. 

We will identify N(p, q) and M(p, q) as the coefficients of certain symmetric functions, 
introduce an operator for extracting coefficients, and show that its action yields the recur- 
sion above. 

Let e n denote the elementary symmetric function of degree n, in a countably infinite num- 
ber of variables {x±, x 2 , • • • }: 

e \ • 7" - /• ... r f 

T\<T2< — <r n 

and let h n be the complete symmetric function of degree n: 

h n . ^ ^ x ri x r2 ■ ■ ■ x rn , 

ri<r2<—<r n 

where r%, ■ ■ ■ , r n G {1, 2, 3, . . . }. For convenience, let x = e = h = 1 and e n = h n = 
if n < 0. Given r G N n , let x r := x^ 1 • • -x^ n and x r := x ri ■ ■ -x Tn . Apply the same 
notation for e r and e r , as well as h r and h r . Note that if r = q, then x r = x q . 
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Lemma 5.3 (MacMahon) For any p e N m , q e N n , 

(1) iV(p, q) w the coefficient ofx p in e q , and 

(2) M(p, q) w the coefficient ofx p in hq. 



Proof The coefficient of x p in e q is the number of ways to choose one term from each of 
the n factors e qi , . . . , e qn , such that the product of these terms is x p . Observe the corre- 
spondence in which the n factors in e q = e qi ■ ■ ■ e qn are identified with the n columns in 
the matrix, and choosing a term x r in a given e qi corresponds to choosing column i to have 
ones in rows r±,...,r qi (and zeros elsewhere). For any choice of terms x r i, . . . , x r ™ from 
e qi , . . . , e qn respectively such that x r i ■ ■ ■ x r n = x p , we have a binary matrix with margins 
(p, q), and conversely, for any such matrix there is such a choice of terms x r i, . . . , x r ™. 
Thus, the coefficient of x p in e q is also the number of such matrices, N(p, q). 

The proof for M(p, q) is the same, except in this case, choosing a term x r in a given 
h q . corresponds to choosing column i to have entries r 1; . . . , r m , and a sequence such that 
x T ■ ■ ■ x r " = x p corresponds to an N- valued matrix with margins (p, q). 



In what follows, when we say "series", we mean a formal power series in x%, x-i y Write 

x = (xi, x 2 , ■ ■ ■) for the sequence of variables, and let Rx = (0, x\, x 2 , ■ ■ ■ ). For k E N, 
define the differential operator: 



i d k 



k\ dx\ 



x=iix 



In other words, after taking the kth derivative with respect to x\ and dividing by k\, replace 

X\ with zero, and with x^ for % — 1, 2, Acting on a series in x x , x 2 , ■ ■ ■ , the operator 

D k annihilates every term except those in which the power of x\ is exactly k. (Note that 
D k coincides with Hammond's operator [fT5l I2T1 . on any symmetric series.) Define 

D T :=D rn .--D n 

(note the reverse order) where n = max{j : rj ^ 0} if r ^ and D r is the identity 
operator otherwise. By applying the operator D r , we keep only terms exactly divisible by 
x T (that is, the power of Xi is for % = 1, . . . , n). In particular, if / is a homogeneous 
series of degree r,, (so that each term has degree ?";), then D r f is a number equal to 
the coefficient of x r in /. Since e q and hq are homogeneous series of degree men by 
Lemma [531 we have 



Corollary 5.4 For any p e N m , q G N n such that YPi = J2 1i> 

(1) D p eq = N(p,q), 

(2) /yi q = M(p,q). | 
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The following identities begin to reveal the utility of the operators D k . 



Lemma 5.5 (MacMahon) For n,k e N, 

(1) D k h, 

(2) D k e r 



(1) D k h n = h n _ k 

e n -k ifk < 1 
if k > 1 
(3) For any functions fi,---,f n > 



sec„(fc) 



Proof (1) and (2) are straightforward calculations. For (3), writing d k = we have 



as? 



k\D k (f 1 ---f n )= £ (d fl h)---(d fn fn) 

re{l,...,n} fe 



x=i?,x 



fc! 



Si! . . . S n ! 



x=Rx 



sGC„(fc) sGC„(fc) 

where the first step follows by recursive application of the product rule, and the second by 
collecting like terms. 



Lemma 5.6 (Power rules) For any k e N, r e N n , 

3 r\s 



sGC r (fc) V 7 

( 2)^= e ( r+ s Ls )^ s - 



sGC* r + is (A;) 

Proof (1) For any m,8 6N, 



by Lemma [531 2 and 3). Thus, 



sec„(fc) 
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(b) r A _n-«i.«i \ r 



E 

sec n (fc) 



Q / r j e n-si+S2 e ^2-S2+s 3 . . , e r„-s„ ^ | r j e r\s 

sec* n (fe) seC r (fc) 

with (a) by Lemma 1531 3). (b) by the preceding observation, (c) by collecting factors, and 
(d) since Q = unless s < r and by the definition of r\s. 

(2) Let m = J2 r i an d let v G N m be any vector such that v = r, so that h r = h v . Let 
S={s:s6 C r+Ls (k)}, and using Lemma l5T27 2), define the surjection g : C v (k) — »■ S by 
g(u) = v' — (v — u)'. Then 

W = J D fc (^ 1 ...^ m )^ ) ^ (AnM • ■ • (A^, J 

uec m (fc) 

= E = E = £ ^ 

ueC ro (fc) ueC v (fc) uec v (fc) 

sec r + £s (fe) ues~ 1 (s) sec r + Ls (fe) ^ ' 

where (a) follows from Lemma I53T 3), (b) by I5.5f l), (c) since hj = if j < and thus 
h v _ u = if u £ v, (d) by|532), and (e) by OP and[52t2). | 

We now complete the generating function proof of Theorem 12.11 If p G N m , q G N n , 
Pi = J2 Qi-> an d r = q, then by Lemma I5T6T 1). 

D p e r = D Lp (D Pl e r ) = ^ f r )D Lp e r \ s , 

seC r ( P i) ^ 

and since e r = e q , then using Corollary 15 .4 1 (twice) we have 

iV(p,r) = iV(p,q) = J D p e q = J D p e r = ^ (A N(Lp, r\s). 

sGC r (pi) ^ ' 

This proves Theorem I2.K 1). Similarly, in view of Corollary 15.41 Theorem I2.1T 2) follows 



immediately from Lemma I5T6T 2). 
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6 Computation time 



Let W(r) := £Li = the weight ofreZ". 

Lemma 6.1 (Properties of the weight) Ifr, s G Z n then 

(1) W(r + s) = W(t) + W(s) 

(2) W(s-Ls) = $> 

(3) W{T\s)=W{r)-J2 8i 

(4) W(s) = $>. 

Proof All four are simple calculations. 

For the rest of this section, fix (p, q) G N m x N n such that J^Pi = Yl ft' an d consider 
(p, q) to be the margins of a set of m x n matrices. First, we address the time to compute 
N(p, q) using Algorithm 12.21 and M(p, q) will follow easily. 

Let V(p, q) denote the set of nontrivial nodes (u, v) in the directed acyclic graph (as dis- 
cussed in Section© descending from (p, q) (including (p, q)), where nontrivial means 
(u,v) 7^ (0,0). Let A k (j) := {s e N k : W(s) = j} for j, k G N. The intuitive content 
of the following lemma is that the graph descending from (p, q) is contained in a union of 
sets A k (j) with weights decreasing by steps of P i, . . . , Pm . 

Lemma 6.2 (Descendants) Z>(p, q) C {(u, v) : u = i^ _1 p, v G A b (tj), j = 1, . . . , m}, 

where tj = Y^iLj Pi an< ^ b = max ft- 

Proof By the form of the recursion, (u, v) G T>(p, q) if and only if for some 1 < j < m 
there exist s 1 , . . . , s 3-1 in C T (pi), ■ ■ ■ , C r3 (Pj-i) respectively, with r 1 = q, r l+1 = r 4 \s* 
for i — 1, . . . , j — 1, such that (u, v) = (L J_1 p, r 3 ). For j > 2, by Lemma [Of 3 and 4), 

W(r*) = H A (i J '- 1 \s i - 1 ) = WV" 1 ) - Pi _i = W^(r^ 2 ) - Pi _ 2 - Pi _i 

n J— 1 m J— 1 

= . . . = ^/(r 1 ) - ( Pl + ■ • • + Pi _ x ) = ^ ft - J^ Pi = J^ Pi - J^ Pi = t j; 

i=l i=l i=l i=l 

and r J G by construction, so r J G A b (tj). Hence, (u, v) = (L J ' _1 p, r 3 ') belongs to the 
set as claimed. 

Let T(p, q) be the time (number of machine operations) required by the algorithm (Algo- 
rithm 12.21) to compute N(p, q) after precomputing all needed binomial coefficients. Let 
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r(u, v) be the time to compute N(vl, v) given N(Lu, v\s) for all s G C v (u\). That is, 
T(p,q) is the time to perform the entire recursive computation, whereas r(u, v) is the 
time to perform a given call to the algorithm not including time spent in subcalls to the 
algorithm. 

Let n := #{i : q { > 0} denote the number of nonempty columns. By constructing 
Pascal's triangle, we precompute all possible binomial coefficients that will be needed, and 
store them in a lookup table. We only need binomial coefficients with entries less or equal 
to no, for the following reason. In the binary case, the recursion involves numbers of the 
form f v ) with s < v, and for any descendent (u, v) and any i > we have Vi < n 
since the number of columns with sum i is less or equal to the total number of nonempty 
columns. For the N- valued case, the same set of binomial coefficients will be sufficient, 
since then we have numbers of the form ( v ^ Ls ) with s < v + Ls, and thus 

Vi + s i+ i < Vi + v i+ i + s i+2 < ■ ■ ■ < Vi + v i+ i + v i+2 H <n , 

where the last inequality holds because the number of columns j with sum greater or equal 
to i is no more than the total number of nonempty columns. Since the addition of two d- 
digit numbers takes 0(d) time, and there are ( n ° 2 f2 ) binomial coefficients with entries less 
or equal to n , then the bound log ({) + 1 < n log 2 + 1 on the number of digits for such a 
binomial coefficient shows that this pre-computation can be done in 0(^q) time. Except in 
trivial cases (when the largest column sum is 1), the additional time needed does not affect 
the bounds on T(p, q) that we will prove below. 

We now bound the time required for a given call to the algorithm. 

Lemma 6.3 (Time per call) r(u,v) < 0((ab + c) (log c) 3 |C f) ('Ui)|) /or (u, v) e T>(p,q), 

where a = maxp^, b = max q it and c = ^2pi. 

Proof Note that we always have < n , since the number of columns with sum i cannot 
exceed the number of nonempty columns. Thus, in the recursion formula, for each s in the 
sum corresponding to (u, v), we have the bound 

Let T m (k) be the time required to multiply two numbers of magnitude k or less. By the 
Schonhage-Strassen algorithm [|30l . T m (k) < 0( (log A;) (log log A;) (log log log k)). There- 
fore, T m (( v )) < T m (c a ) < 0(a(logc) 3 ). Since we have precomputed the binomial co- 
efficients, the time required to compute ( v ) is thus bounded by 0(a6(logc) 3 ). To finish 
computing the term corresponding to s in the recursion formula, we must multiply ( v ) by 
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N(Lu, v\s). Since 



iV(Lu,v\s) < iV(p,q) < ] [ ( n ° ) < J]< = n c < c c , 

i=i VP*/ i=i 

then this multiplication can be done in T m (N(p, q)) < T m (c c ) < 0(c(logc) 3 ) time. Since 
we are summing over C v (ui), and C v (ui) C Cb{u\), then altogether we have r(u, v) < 
0((afe + c)(logc) 3 |Cb(Mi)|) for the time per call. 

(j + k — 1\ 
for any j, fcGN. 
fc-i y 

Proof The map /(r) = (lri, 2r 2 , . . . , krk) is an injection / : A k (j) — >■ C k (j). Thus, 

#A fc ( J )<#c fc (j) = e't!7 1 ). I 



We are now ready to prove Theorem 12. 3 1 for the case of N(p, q). 



Proof of Theorem H3] for N(p, q) 

By storing intermediate results in a lookup table, once we have computed _/V(u, v) upon 
our first visit to node (u, v), we can simply reuse the result for later visits. Hence, we need 
only expend r(u, v) time for each node (u, v) occuring in the graph. Let tj = Y^iLj Pi an d 
d = (ab + c)(logc) 3 . Then 

(a) m 

T( P ,q)= Yl ^( u ^)<E E < Lj ~^) 

(u,v)ex>(p,q) i=i veAsfe) 

< EE ^!^)!) = E°(^fe)ll A ^;)i) 

j v j 



<E°( d 



p i + 6-l\ /t i + 6-l 
b-1 A 6-1 



6-1^ 



where (a) follows by Lemma |6\21 (b) by 16.31 (c) by 16.41 and (d) since pj < a and tj < c. 
This proves (1) and (2). Now, (3) and (4) follow from (2) since a < c < bn. 
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Proof of Theorem |23] for M(p, q) 



Other than the coefficients, the only difference between the recursion for M(p, q) and 
that for N(p, q) is that we are summing over s such that s G C T+Ls (pi). Lemma I6T21 
holds with the same proof, except with C T (pi), ■ ■ ■ , C rl (Pj~i) replaced by C r +Ls (pi), 
. . . , C rJ 1+isJ (pj-i), respectively. Considering Lemma l63l let (u, v) be a descendent 
of (p, q) in the graph for M(p, q), and let s be such that s G C v+Ls (ui). Similarly to 
before, recalling that v t + s i+1 < n (as proven above in our discussion of precomputing 
the binomial coefficients), we have 

C t Ls ) = n + s f +1 ) * n<«. + s »o E " < < < 

This yields T m (f + s Ls )) < T m {c a ) < 0(a(logc) 3 ), just as before. Since 

m(lu, v\ S ) < m( P , q ) < n ( p ° j < ri( 2c ) pi = ( 2c ) C ' 

i=l ^ ' j 

then we also obtain T m (M(p, q)) < T m ((2c) c ) < 0(c(logc) 3 ) as before. Further, {s : s G 
C v+Ls (mi)} C C b {ui), so altogether the time per call is 0{{ab + c)(logc) 3 |C 6 (Mi)|), and 
thus the result of Lemma [631 continues to hold. With this result, the proof of the bounds 
goes through as well. 

This completes the proof of Theorem 12.31 Now, we address the time required to uniformly 
sample a matrix with specified margins. Let T r (k) be the maximum over 1 < j < k 
of the expected time to generate a random integer uniformly between 1 and j. If we are 
given a random bitstream (independent and identically distributed Bernoulli(l/2) random 
variables) with constant cost per bit, then T r ( k) = O ( log k), since for any j < k, |~log 2 j] < 
[log 2 k~\ random bits can be used to generate an integer uniformly between 1 and 2^ og2 ^ 
and then rejection sampling can be used to generate uniform samples over {1, . . . , j}. Since 
the expected value of a Geometric (p) random variable is 1/p, then the expected number of 
samples required to obtain one that falls in {1, . . . , j} is always less than 2. More generally, 
for any fixed d G N, if we can draw uniform samples from {1, . . . , d}, then we have 
T r {k) = 0(log k) by considering the base-rf analogue of the preceding argument. 

Lemma 6.5 (Sampling time) A Igorithm \2.4\ takes 0(mT r (n c ) +maT r (n) +mb log(a + 6)) 
expected time per sample in the binary case, and 0(mT r ((2c) c ) + maT r {n) + mMog(a + 
b)) expected time per sample in the N-valued case. If T r (k) = O(logfc), then this is 
0(mc log c) expected time per sample in both cases. 

Remark If b is bounded then O (mc log c) < O(mnlogn) since c < bn, and so this is 
polynomial expected time for bounded column sums. 
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Proof By the form of the recursion, the depth of the graph descending from (p, q) is equal 
to the number of rows m, since p G N m and thus L m p = 0. For each of the m iterations 
of the sampling algorithm, we begin at some node (u, v), and we must (A) randomly 
choose a child (Lu, v\s) with probability proportional to its count times the number of 
corresponding rows, and then (B) choose a row uniformly from among the Q) possible 
choices in the binary case (or ( v ^ is ) in the N- valued case). 

First consider the binary case. To randomly choose a child, consider a partition of the 
integers 1, . . . , N(u, v) with each part corresponding to a term in the recursion formula 
for N(u, v). Generate an integer uniformly at random between 1 and N(u, v), and 
choose the corresponding child. Generating such a random number takes T r (N(u, v)) < 
T r (iV(p,q)) < T r (n c ) time. Since there are no more than ( a ^ 1 ) < (a + 6 - 
children at any step, one can determine which child corresponds to the chosen number in 
0((6 — 1) log(a + 6 — 1)) time by organizing the children in a binary tree. So (A) takes 
0(T r (riQ) + 61og(a + 6)) time. Choosing a row consists of uniformly sampling a subset of 
size Si from a set of Vi elements, for i — 1, . . . , 6. Sampling such a subset can be done by 
sampling without replacement Sj times, which takes YljLo T r (vi — j) < SiT r (n ) time. So 

(B) can be done in JZ i=1 SjT r (n ) < aT r (n ) time. Repeating this process m times, once 
for each row, we see that sampling a matrix takes 0{mT r {n^) +m61og(a + 6) + maT r (n )) 
time. IfT r (k) < 0(log k), this is 0(mc log n + m61og(a + 6) + ma log n ) < O(mclogc) 
since a, 6, n < c. 

For the N-valued case, the same argument applies, replacing N(u, v) with M(u, v), Uq 
with (2c) c , and V{ with + s i+ x- 
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A Enumeration results 

Binary matrices with margins (70, 30, 20, 10, 5®, 3^,2^), (4( 80 \ 3< 20 >) 

860585058801817078819959949756041558231879514104670757612387 
280341919502865086909993523205599348663646837362726765460951 
032776118129432733489342067673016169716787054236343091407458 
802261593735765113169808512677339861494709092492858489355535 
514748397544147637928475318462070009855280569561693514768239 
201499080842592443823774161366680107327323365049702068246736 
456919918589686056321467354298509024976141650428747522863473 
5295 1 52693 1 8246400000000000000000000000 



N-valued matrices with margins (70, 30, 20, 10, 5< 6 \ 4< 10 ), 3^ 20 \ 2< 60 )), (4( 80 ), 3( 20 )) 

620017488391049592297896956531192562528805388295441812965295 
130897484012791595142882674755488640101825726867156331426482 
441148514978852842582445295040041143220637964258279947442682 
896809706562683189375098411751981435132377208717294759756041 
358372207736032818841045369779439398975681041714752821787419 
816573563436066161167632677774184809010338787868042742993719 
703936093873250600121874335524794990013547042810153560084573 
1 3303573 12 17642637607 1 536 1 56 1 1 02985 1 392000000000000000000000 
000 



Ehrhart polynomials H n (r) for n — 4, . . . , 8 

H 4 (l) = 24 
# 4 (2) = 282 
if 4 (3) = 2008 

H 5 (l) = 120 
H 5 (2) = 6210 
# 5 (3) = 153040 
H b {A) = 2224955 
#5(5) = 22069251 
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H 5 (6) = 164176640 



H 6 


,1) = 


720 


TT 

H 6 


2) = 


202410 


Hp. 




20933840 


H 6 


;4) = 


1047649905 


H Q 


;5) = 


30767936616 


H 6 


;e) = 


602351808741 


H Q 


7) = 


8575979362560 


H e 


:§) = 


94459713879600 


H 6 


[9) = 


842286559093240 


H 6 


;io) = 


= 6292583664553881 


H 7 


,i) = 


5040 


TT 

H 7 


2) = 


9135630 


TT 

H 7 


3) = 


4662857360 


TT 

H 7 


' A \ 

4) = 


936670590450 


TT 

H 7 


5) = 


94161778046406 


TT 

H 7 


6) = 


5562418293759978 


TT 

H 7 


J) = 


21571760804651 1873 


n 7 


'%) - 


5945968652327831 925 


H 7 


;s) = 


123538613356253145400 


H r 


10) = 


= 2023270039486328373811 


H 7 


;n) = 


= 27046306550096288483238 


H 7 


;i2) = 


= 303378141987182515342992 


H 7 


;i3) = 


= 2920054336492521720572276 


H 7 


;i4) = 


= 24563127009195223721952590 


H 7 


;i5) = 


= 183343273080700916973016745 


H 8 


' -i \ 

,i) = 


40320 


TT 

H 8 


2) = 


545007960 


TT 

H 8 


3) = 


1579060246400 


TT 

H 8 


f A \ 

4) = 


1455918295922650 


TT 

H 8 


,5) = 


569304690994400256 


TT 

H 8 


6) = 


H4601242382721619224 


TT 

H 8 


J) = 


loo9U 4ly4zo4zzo4o9U4 


H 8 


[8) = 


1046591482728407939338275 


H 8 


[9) = 


56272722406349235035916800 


H 8 


;io) = 


= 2233160342369825596702148720 


H 8 


;n) = 


= 68316292103293669997188919040 


H 8 


;i2) = 


= 1667932098862773837734823042196 


H 8 


13), 


= 33427469280977307618866364694400 


H 8 


:i4) = 


= 562798805673342016752366344185200 
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H 8 (15) = 8115208977465404874100226492575360 
H 8 (16) = 101857066150530294146428615917957029 
H 8 (17) = 1128282526405022554049557329097252992 
H 8 (18) = 11161302946841260178530673680176000200 
#8(19) = 99613494890126594335550124219924540800 
H 8 (20) = 809256770610540675454657517194018680846 
H 8 (21) = 6031107989875562751266116901999327710720 
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